MCMC convergence diagnostics

Abstract

Background In the report naomi-simple_fit with parameter tmbstan = TRUE, we used the NUTS algorithm to perform MCMC inference for the simplified Naomi model.

Task Here we assess whether or not the results of the MCMC are suitable using a range of diagnostic tools.

We start by obtaining results from the latest version of naomi-simple_fit with tmbstan = TRUE.

out <- readRDS("depends/out.rds")
mcmc <- out$mcmc$stanfit

This MCMC took 3.28033377378351 to run

bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())

1 \(\hat R\)

We are looking for values of \(\hat R\) less than 1.1 here.

rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat(rhats)

(big_rhats <- rhats[rhats > 1.1])

## named numeric(0)

length(big_rhats) / length(rhats)

## [1] 0

2 ESS ratio

Reasonable to be worried about values less than 0.1 here.

ratios <- bayesplot::neff_ratio(mcmc)
bayesplot::mcmc_neff(ratios)

3 ESS

What are the total effective sample sizes?

#' I think that this $summary should be all of the chains grouped together
mcmc_summary <- summary(mcmc)$summary

data.frame(mcmc_summary) %>%
  tibble::rownames_to_column("param") %>%
  ggplot(aes(x = n_eff)) +
    geom_histogram(alpha = 0.8) +
    labs(x = "ESS", y = "Count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

4 Autocorrelation

How much autocorrelation is there in the chains?

bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))

5 Univariate traceplots

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))

5.1 Prevalence model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))

5.2 ART model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))

5.3 Other

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))

5.4 ANC model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model

6 Pairs plots

There is a prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable. Let’s have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter.

area_merged <- sf::read_sf(system.file("extdata/demo_areas.geojson", package = "naomi"))

nb <- area_merged %>%
  filter(area_level == max(area_level)) %>%
  bsae::sf_to_nb()

neighbours_pairs_plot <- function(par, i) {
  neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
  bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}

# area_merged %>%
#   filter(area_level == max(area_level)) %>%
#   print(n = Inf)

Here are Nkhata Bay and neighbours:

neighbours_pairs_plot("log_or_gamma", 5) 

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## ℹ The deprecated feature was likely used in the bayesplot package.
##   Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

And here are Blantyre and neighbours:

neighbours_pairs_plot("log_or_gamma", 26)

7 NUTS specific assessment

np <- bayesplot::nuts_params(mcmc)

Are there any divergent transitions?

np %>%
  filter(Parameter == "divergent__") %>%
  summarise(n_divergent = sum(Value))

bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))

We can also use energy plots (Betancourt 2017): ideally these two histograms would be the same When the histograms are quite different, it may suggest the chains are not fully exploring the tails of the target distribution.

bayesplot::mcmc_nuts_energy(np)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Betancourt, Michael. 2017. “A Conceptual Introduction to Hamiltonian Monte Carlo.” arXiv Preprint arXiv:1701.02434.